home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dtrsen.f < prev    next >
Text File  |  1996-07-19  |  14KB  |  422 lines
  1.       SUBROUTINE DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
  2.      $                   M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          COMPQ, JOB
  11.       INTEGER            INFO, LDQ, LDT, LIWORK, LWORK, M, N
  12.       DOUBLE PRECISION   S, SEP
  13. *     ..
  14. *     .. Array Arguments ..
  15.       LOGICAL            SELECT( * )
  16.       INTEGER            IWORK( * )
  17.       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ),
  18.      $                   WR( * )
  19. *     ..
  20. *
  21. *  Purpose
  22. *  =======
  23. *
  24. *  DTRSEN reorders the real Schur factorization of a real matrix
  25. *  A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in
  26. *  the leading diagonal blocks of the upper quasi-triangular matrix T,
  27. *  and the leading columns of Q form an orthonormal basis of the
  28. *  corresponding right invariant subspace.
  29. *
  30. *  Optionally the routine computes the reciprocal condition numbers of
  31. *  the cluster of eigenvalues and/or the invariant subspace.
  32. *
  33. *  T must be in Schur canonical form (as returned by DHSEQR), that is,
  34. *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
  35. *  2-by-2 diagonal block has its diagonal elemnts equal and its
  36. *  off-diagonal elements of opposite sign.
  37. *
  38. *  Arguments
  39. *  =========
  40. *
  41. *  JOB     (input) CHARACTER*1
  42. *          Specifies whether condition numbers are required for the
  43. *          cluster of eigenvalues (S) or the invariant subspace (SEP):
  44. *          = 'N': none;
  45. *          = 'E': for eigenvalues only (S);
  46. *          = 'V': for invariant subspace only (SEP);
  47. *          = 'B': for both eigenvalues and invariant subspace (S and
  48. *                 SEP).
  49. *
  50. *  COMPQ   (input) CHARACTER*1
  51. *          = 'V': update the matrix Q of Schur vectors;
  52. *          = 'N': do not update Q.
  53. *
  54. *  SELECT  (input) LOGICAL array, dimension (N)
  55. *          SELECT specifies the eigenvalues in the selected cluster. To
  56. *          select a real eigenvalue w(j), SELECT(j) must be set to
  57. *          .TRUE.. To select a complex conjugate pair of eigenvalues
  58. *          w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
  59. *          either SELECT(j) or SELECT(j+1) or both must be set to
  60. *          .TRUE.; a complex conjugate pair of eigenvalues must be
  61. *          either both included in the cluster or both excluded.
  62. *
  63. *  N       (input) INTEGER
  64. *          The order of the matrix T. N >= 0.
  65. *
  66. *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
  67. *          On entry, the upper quasi-triangular matrix T, in Schur
  68. *          canonical form.
  69. *          On exit, T is overwritten by the reordered matrix T, again in
  70. *          Schur canonical form, with the selected eigenvalues in the
  71. *          leading diagonal blocks.
  72. *
  73. *  LDT     (input) INTEGER
  74. *          The leading dimension of the array T. LDT >= max(1,N).
  75. *
  76. *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
  77. *          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
  78. *          On exit, if COMPQ = 'V', Q has been postmultiplied by the
  79. *          orthogonal transformation matrix which reorders T; the
  80. *          leading M columns of Q form an orthonormal basis for the
  81. *          specified invariant subspace.
  82. *          If COMPQ = 'N', Q is not referenced.
  83. *
  84. *  LDQ     (input) INTEGER
  85. *          The leading dimension of the array Q.
  86. *          LDQ >= 1; and if COMPQ = 'V', LDQ >= N.
  87. *
  88. *  WR      (output) DOUBLE PRECISION array, dimension (N)
  89. *  WI      (output) DOUBLE PRECISION array, dimension (N)
  90. *          The real and imaginary parts, respectively, of the reordered
  91. *          eigenvalues of T. The eigenvalues are stored in the same
  92. *          order as on the diagonal of T, with WR(i) = T(i,i) and, if
  93. *          T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and
  94. *          WI(i+1) = -WI(i). Note that if a complex eigenvalue is
  95. *          sufficiently ill-conditioned, then its value may differ
  96. *          significantly from its value before reordering.
  97. *
  98. *  M       (output) INTEGER
  99. *          The dimension of the specified invariant subspace.
  100. *          0 < = M <= N.
  101. *
  102. *  S       (output) DOUBLE PRECISION
  103. *          If JOB = 'E' or 'B', S is a lower bound on the reciprocal
  104. *          condition number for the selected cluster of eigenvalues.
  105. *          S cannot underestimate the true reciprocal condition number
  106. *          by more than a factor of sqrt(N). If M = 0 or N, S = 1.
  107. *          If JOB = 'N' or 'V', S is not referenced.
  108. *
  109. *  SEP     (output) DOUBLE PRECISION
  110. *          If JOB = 'V' or 'B', SEP is the estimated reciprocal
  111. *          condition number of the specified invariant subspace. If
  112. *          M = 0 or N, SEP = norm(T).
  113. *          If JOB = 'N' or 'E', SEP is not referenced.
  114. *
  115. *  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
  116. *
  117. *  LWORK   (input) INTEGER
  118. *          The dimension of the array WORK.
  119. *          If JOB = 'N', LWORK >= max(1,N);
  120. *          if JOB = 'E', LWORK >= M*(N-M);
  121. *          if JOB = 'V' or 'B', LWORK >= 2*M*(N-M).
  122. *
  123. *  IWORK   (workspace) INTEGER array, dimension (LIWORK)
  124. *          IF JOB = 'N' or 'E', IWORK is not referenced.
  125. *
  126. *  LIWORK  (input) INTEGER
  127. *          The dimension of the array IWORK.
  128. *          If JOB = 'N' or 'E', LIWORK >= 1;
  129. *          if JOB = 'V' or 'B', LIWORK >= M*(N-M).
  130. *
  131. *  INFO    (output) INTEGER
  132. *          = 0: successful exit
  133. *          < 0: if INFO = -i, the i-th argument had an illegal value
  134. *          = 1: reordering of T failed because some eigenvalues are too
  135. *               close to separate (the problem is very ill-conditioned);
  136. *               T may have been partially reordered, and WR and WI
  137. *               contain the eigenvalues in the same order as in T; S and
  138. *               SEP (if requested) are set to zero.
  139. *
  140. *  Further Details
  141. *  ===============
  142. *
  143. *  DTRSEN first collects the selected eigenvalues by computing an
  144. *  orthogonal transformation Z to move them to the top left corner of T.
  145. *  In other words, the selected eigenvalues are the eigenvalues of T11
  146. *  in:
  147. *
  148. *                Z'*T*Z = ( T11 T12 ) n1
  149. *                         (  0  T22 ) n2
  150. *                            n1  n2
  151. *
  152. *  where N = n1+n2 and Z' means the transpose of Z. The first n1 columns
  153. *  of Z span the specified invariant subspace of T.
  154. *
  155. *  If T has been obtained from the real Schur factorization of a matrix
  156. *  A = Q*T*Q', then the reordered real Schur factorization of A is given
  157. *  by A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span
  158. *  the corresponding invariant subspace of A.
  159. *
  160. *  The reciprocal condition number of the average of the eigenvalues of
  161. *  T11 may be returned in S. S lies between 0 (very badly conditioned)
  162. *  and 1 (very well conditioned). It is computed as follows. First we
  163. *  compute R so that
  164. *
  165. *                         P = ( I  R ) n1
  166. *                             ( 0  0 ) n2
  167. *                               n1 n2
  168. *
  169. *  is the projector on the invariant subspace associated with T11.
  170. *  R is the solution of the Sylvester equation:
  171. *
  172. *                        T11*R - R*T22 = T12.
  173. *
  174. *  Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote
  175. *  the two-norm of M. Then S is computed as the lower bound
  176. *
  177. *                      (1 + F-norm(R)**2)**(-1/2)
  178. *
  179. *  on the reciprocal of 2-norm(P), the true reciprocal condition number.
  180. *  S cannot underestimate 1 / 2-norm(P) by more than a factor of
  181. *  sqrt(N).
  182. *
  183. *  An approximate error bound for the computed average of the
  184. *  eigenvalues of T11 is
  185. *
  186. *                         EPS * norm(T) / S
  187. *
  188. *  where EPS is the machine precision.
  189. *
  190. *  The reciprocal condition number of the right invariant subspace
  191. *  spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP.
  192. *  SEP is defined as the separation of T11 and T22:
  193. *
  194. *                     sep( T11, T22 ) = sigma-min( C )
  195. *
  196. *  where sigma-min(C) is the smallest singular value of the
  197. *  n1*n2-by-n1*n2 matrix
  198. *
  199. *     C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
  200. *
  201. *  I(m) is an m by m identity matrix, and kprod denotes the Kronecker
  202. *  product. We estimate sigma-min(C) by the reciprocal of an estimate of
  203. *  the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C)
  204. *  cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
  205. *
  206. *  When SEP is small, small changes in T can cause large changes in
  207. *  the invariant subspace. An approximate bound on the maximum angular
  208. *  error in the computed right invariant subspace is
  209. *
  210. *                      EPS * norm(T) / SEP
  211. *
  212. *  =====================================================================
  213. *
  214. *     .. Parameters ..
  215.       DOUBLE PRECISION   ZERO, ONE
  216.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  217. *     ..
  218. *     .. Local Scalars ..
  219.       LOGICAL            PAIR, SWAP, WANTBH, WANTQ, WANTS, WANTSP
  220.       INTEGER            IERR, K, KASE, KK, KS, N1, N2, NN
  221.       DOUBLE PRECISION   EST, RNORM, SCALE
  222. *     ..
  223. *     .. External Functions ..
  224.       LOGICAL            LSAME
  225.       DOUBLE PRECISION   DLANGE
  226.       EXTERNAL           LSAME, DLANGE
  227. *     ..
  228. *     .. External Subroutines ..
  229.       EXTERNAL           DLACON, DLACPY, DTREXC, DTRSYL, XERBLA
  230. *     ..
  231. *     .. Intrinsic Functions ..
  232.       INTRINSIC          ABS, MAX, SQRT
  233. *     ..
  234. *     .. Executable Statements ..
  235. *
  236. *     Decode and test the input parameters
  237. *
  238.       WANTBH = LSAME( JOB, 'B' )
  239.       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
  240.       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
  241.       WANTQ = LSAME( COMPQ, 'V' )
  242. *
  243.       INFO = 0
  244.       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.WANTS .AND. .NOT.WANTSP )
  245.      $     THEN
  246.          INFO = -1
  247.       ELSE IF( .NOT.LSAME( COMPQ, 'N' ) .AND. .NOT.WANTQ ) THEN
  248.          INFO = -2
  249.       ELSE IF( N.LT.0 ) THEN
  250.          INFO = -4
  251.       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  252.          INFO = -6
  253.       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
  254.          INFO = -8
  255.       ELSE
  256. *
  257. *        Set M to the dimension of the specified invariant subspace,
  258. *        and test LWORK and LIWORK.
  259. *
  260.          M = 0
  261.          PAIR = .FALSE.
  262.          DO 10 K = 1, N
  263.             IF( PAIR ) THEN
  264.                PAIR = .FALSE.
  265.             ELSE
  266.                IF( K.LT.N ) THEN
  267.                   IF( T( K+1, K ).EQ.ZERO ) THEN
  268.                      IF( SELECT( K ) )
  269.      $                  M = M + 1
  270.                   ELSE
  271.                      PAIR = .TRUE.
  272.                      IF( SELECT( K ) .OR. SELECT( K+1 ) )
  273.      $                  M = M + 2
  274.                   END IF
  275.                ELSE
  276.                   IF( SELECT( N ) )
  277.      $               M = M + 1
  278.                END IF
  279.             END IF
  280.    10    CONTINUE
  281. *
  282.          N1 = M
  283.          N2 = N - M
  284.          NN = N1*N2
  285. *
  286.          IF( LWORK.LT.1 .OR. ( ( WANTS .AND. .NOT.WANTSP ) .AND.
  287.      $       LWORK.LT.NN ) .OR. ( WANTSP .AND. LWORK.LT.2*NN ) ) THEN
  288.             INFO = -15
  289.          ELSE IF( LIWORK.LT.1 .OR. ( WANTSP .AND. LIWORK.LT.NN ) ) THEN
  290.             INFO = -17
  291.          END IF
  292.       END IF
  293.       IF( INFO.NE.0 ) THEN
  294.          CALL XERBLA( 'DTRSEN', -INFO )
  295.          RETURN
  296.       END IF
  297. *
  298. *     Quick return if possible.
  299. *
  300.       IF( M.EQ.N .OR. M.EQ.0 ) THEN
  301.          IF( WANTS )
  302.      $      S = ONE
  303.          IF( WANTSP )
  304.      $      SEP = DLANGE( '1', N, N, T, LDT, WORK )
  305.          GO TO 40
  306.       END IF
  307. *
  308. *     Collect the selected blocks at the top-left corner of T.
  309. *
  310.       KS = 0
  311.       PAIR = .FALSE.
  312.       DO 20 K = 1, N
  313.          IF( PAIR ) THEN
  314.             PAIR = .FALSE.
  315.          ELSE
  316.             SWAP = SELECT( K )
  317.             IF( K.LT.N ) THEN
  318.                IF( T( K+1, K ).NE.ZERO ) THEN
  319.                   PAIR = .TRUE.
  320.                   SWAP = SWAP .OR. SELECT( K+1 )
  321.                END IF
  322.             END IF
  323.             IF( SWAP ) THEN
  324.                KS = KS + 1
  325. *
  326. *              Swap the K-th block to position KS.
  327. *
  328.                IERR = 0
  329.                KK = K
  330.                IF( K.NE.KS )
  331.      $            CALL DTREXC( COMPQ, N, T, LDT, Q, LDQ, KK, KS, WORK,
  332.      $                         IERR )
  333.                IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
  334. *
  335. *                 Blocks too close to swap: exit.
  336. *
  337.                   INFO = 1
  338.                   IF( WANTS )
  339.      $               S = ZERO
  340.                   IF( WANTSP )
  341.      $               SEP = ZERO
  342.                   GO TO 40
  343.                END IF
  344.                IF( PAIR )
  345.      $            KS = KS + 1
  346.             END IF
  347.          END IF
  348.    20 CONTINUE
  349. *
  350.       IF( WANTS ) THEN
  351. *
  352. *        Solve Sylvester equation for R:
  353. *
  354. *           T11*R - R*T22 = scale*T12
  355. *
  356.          CALL DLACPY( 'F', N1, N2, T( 1, N1+1 ), LDT, WORK, N1 )
  357.          CALL DTRSYL( 'N', 'N', -1, N1, N2, T, LDT, T( N1+1, N1+1 ),
  358.      $                LDT, WORK, N1, SCALE, IERR )
  359. *
  360. *        Estimate the reciprocal of the condition number of the cluster
  361. *        of eigenvalues.
  362. *
  363.          RNORM = DLANGE( 'F', N1, N2, WORK, N1, WORK )
  364.          IF( RNORM.EQ.ZERO ) THEN
  365.             S = ONE
  366.          ELSE
  367.             S = SCALE / ( SQRT( SCALE*SCALE / RNORM+RNORM )*
  368.      $          SQRT( RNORM ) )
  369.          END IF
  370.       END IF
  371. *
  372.       IF( WANTSP ) THEN
  373. *
  374. *        Estimate sep(T11,T22).
  375. *
  376.          EST = ZERO
  377.          KASE = 0
  378.    30    CONTINUE
  379.          CALL DLACON( NN, WORK( NN+1 ), WORK, IWORK, EST, KASE )
  380.          IF( KASE.NE.0 ) THEN
  381.             IF( KASE.EQ.1 ) THEN
  382. *
  383. *              Solve  T11*R - R*T22 = scale*X.
  384. *
  385.                CALL DTRSYL( 'N', 'N', -1, N1, N2, T, LDT,
  386.      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
  387.      $                      IERR )
  388.             ELSE
  389. *
  390. *              Solve  T11'*R - R*T22' = scale*X.
  391. *
  392.                CALL DTRSYL( 'T', 'T', -1, N1, N2, T, LDT,
  393.      $                      T( N1+1, N1+1 ), LDT, WORK, N1, SCALE,
  394.      $                      IERR )
  395.             END IF
  396.             GO TO 30
  397.          END IF
  398. *
  399.          SEP = SCALE / EST
  400.       END IF
  401. *
  402.    40 CONTINUE
  403. *
  404. *     Store the output eigenvalues in WR and WI.
  405. *
  406.       DO 50 K = 1, N
  407.          WR( K ) = T( K, K )
  408.          WI( K ) = ZERO
  409.    50 CONTINUE
  410.       DO 60 K = 1, N - 1
  411.          IF( T( K+1, K ).NE.ZERO ) THEN
  412.             WI( K ) = SQRT( ABS( T( K, K+1 ) ) )*
  413.      $                SQRT( ABS( T( K+1, K ) ) )
  414.             WI( K+1 ) = -WI( K )
  415.          END IF
  416.    60 CONTINUE
  417.       RETURN
  418. *
  419. *     End of DTRSEN
  420. *
  421.       END
  422.